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FAST ALGORITHMS FOR COMBUSTION KINETICS CALCULATIONS: A COMPARISON* 

Krishnan Radhakrishnan** 

NASA Lewis Research Center 


Many practical problems arising in chemically reacting flows require the 
simultaneous numerical integration of large sets of chemical kinetic rate 
equations of the type shown in figure 1. The initial value problem is that of 
finding the composition and temperature at the end of a prescribed time 
interval, given the initial mixture composition and temperature, the pressure, 
and the reaction mechanism. Multi-dimensional modeling of reactive flows 
requires the integration of the system of ordinary differential equations 
(ode's) given in figure 1 at several thousand grid points. To make such 
calculations practicable, it is necessary to have a very fast batch chemistry 
integrator . 

T" identify the fastest algorithm currently available for the numerical 
integration of chemical kinetic rat ? equations, several algorithms have been 
examined. In the present paper, we summarize our findings to date — details 
are available in references (1) and (2). The algorithms examined in this work 
include t*a general-purpose codes EPISODE and LSODE (refs. 3 and 4), and three 
special-purpose (for chemical kinetic calculations) codes CHEMEQ (ref. 5), 
CREK1D (refs. 6 and 7), and GCKP84 (refs. 8 and 9). In addition, an explicit 
Runge-Kutta-Merson differential equation solver (ref. 10) (IMSL Routine DASCRU) 
is used to illustrate the problems associated with integrating che - "ml kinetic 
rate equations by a classical method. These methods are summarized in figure 3. 

The algorithms summarized in figure 3 were applied to two test problems 
drawn from combustion kinetics. These problems, summarized in figure 4, 
included all three combustion regimes: induction, heat release ar ‘ 
equilibration. Figures b and 6 present variations of the temperature and 
species mole fractions with time for test problems 1 and 2, respectively. Both 
test problems were integrated over a time internal c' 1 1 ms in order to obtain 
near-equilibration of all species and temperature. 

Of the codes examined in this study, only CUEK1D and GCKP84 were written 
explicitly for integrating exothermic, non- isothermal combustion rate 
equations. These therefore have built-in procedures for calculating the 
temperature (T). For the other codes, two different methods, labeled as Methods 
A and B, were used to compute T. The following convention was adopted in naming 
these other codes: those using temperature method A were given the suffix-A 
(e.g. LSODE-A) and those using temperature method B were given the suffix-B 
(e.g.. LSODE-B). In Method A, T was calculated from the mole numbers and the 
initial mixture enthalpy using an algebraic energy conservation equation (given 
in figure 7) and a Newton-Raphson iteration technique. In this method, the 
temperature is not an explicit independent variable, so the number of ode's is 
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equal to the number, NS, of distinct chemical species in the mixture. The 
integrator therefore tracks only the solutions for the species mole numbers. In 
Method B, the temperature was treated as an additional independent variable and 
evaluated by integrating its time-derivative given in figure 7. In this method, 
the number of ode’s is equal to NS+1, and the integrator tracks the soluti ’ s 
for both the temperature and the species mole numbers. 

All codes were run on the NASA Lewis Research Center’s IBM 370/3033 
computer using single-precision accuracy, except GCKP84 which was in double- 
precision. A typical computational run consisted of initializing the species 
mole numbers, temperature, and CPU time. The integrator was then called with 
values for the necessary input parameters. On return from the integrator, the 
total CPU time required to solve the test problem was calculated. Other 
performance parameters were also recorded — see reference (2) for details. 

Figures 8 and 9 present the computational work (expressed as the CPU time 
in seconds) plotted against the local error tolerance, EPS, for test problems 1 
and 2, respectively. For all codes except EPISODE, EPS is the local relative 
error tolerance. For EPISODE, EPS is a mixed error tolerance — relative fo^ 
species with initially nonzero mole numbers and for the temperature (method B) 
and absolute for species with initially zero mole numbers. Also shown on 
figures 8 and 9 are the CPU times required by the explicit Runge-Kutta method 
for one value of EPS. Note the excessive CPU times required by this technique. 
Its use would make multidimensional modeling of practical combustion devices 
prohibitively expensive. 

For test problem 1, very small values for EPS had to be used for EPISODE 
(figure 8). For values of EPS i 5 x 10“*, EPISODE predicted little or no change 
in the composition and temperature after an elapsed time of 1 ms. Similar 
remarks apply to test problem 2 (figure 9), for which values of 10“'* and 10“ 3 
had to be used for EPISODE-A and EPISODE-B, respectively. Although the runs 
with EPISODE-B and EPS £ 5 x 10““ were successfully completed, the solutions 
(especially for minor species) were significantly different from those given in 
figure 6. With GCKP84 and EPS = 10“ 2 , the solution for test problem 1 exhibited 
serious instability and so was terminated. A more detailed discussion of the 
accuracy of the codes tested in this study can be found in reference (2). 

Examination of figure 8 shows that the difference in computational work 
required by methods A and B is small for test problem 1, with method B being 
more efficient. For test problem 2 (figure 9), the difference is small for 
large values of EPS. But for small values of EPS the difference is more marked, 
with method A being signiflcanlty superior to method B. 

Figures 8 and 9 show that LSODE and CREK1D are superior to the other 
codes. EPISODE is an attractive alternative, especially for test probln 2. 
However, in using EPISODE, a word of caution is in order. The computational 
work can bo strongly dependent on the value for the initial steplength (HO) 
•elected by the user. A poor guess for HO can make EPISODE prohibitively 
expensive to use. Figure 10 Illustrates this behavior for test prvMem 2. Note 
an order of magnitude Increase in the CPU time for a change in HO from 10" 7 to 
10 a * s. Although not shown here, a poor guess for HO also resulted in 
Inaccurate and unstable solutions. In addition, as discussed in refeienee (2), 
the error control performed by EPISODE is unsatisfactory for problems of the 



type eyamined in this study. 

A simple method for increasing the efficiency of the algorithms as applied 
to the present problem was explored. This involved updating the rate constants 
kj and k.j (which was calculated from kj and the concentration equilibrium 
constant) only for temperature changes greater than an amount AT. To avoid a 
trial and error search for the optimum value of AT — defined as that value 
which results in minimum computational work — an approximation for it was 
derived and is presented in figure 11. Comparisons of figure 12 with figures 7 
and 8 show the significant reductions in computational work realized by use of 
the above approximation for AT. 
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OF POOR QUALITY 

Adiabatic, Constant-Pressure, Gas-Phase Chemical Reaction 
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molar rate of formation of species i per unit .lass of mixture, 
kmol-i/kg-mixture s 

forward and reverse rate constants for reaction j 
mole number of species i, kmol-i/kg-mixture 

pre-exponential constants in forward and reverse rate equations for 
reaction j 

activation energy in forward and reverse rate equations for reaction j, 
cal /mole 

number of distinct elementary reactions in mechanism 
number of distinct species in gas mixture 
universal gas constant, 1.987 cal/mol K 

forward and reverse molar reaction rates per unit volume for reaction j, 

kmol/m 3 s 

temperature, K 

mixture mass-density, kg/m 3 

stoichiometric coefficients of species i in reaction j as a reactant, 
and as a product, respectively 


Figure 1 Governing Ordinary Differential Equations 


260 




: • t 

i f 


Given 

Initial Mixture Composition and Ter.^rature 
Pressure 

Reaction Mechanism 

Find, at the End of a Prescribed Time Interval 
Mixture Composition and Temperature 

Figure 2 Problem Statement 


Method 

Description 

CCKP84 

Details not yet available. 

r.EKlD 

Variable-step, predictor-corrector method based on an 
exponentially-fitted trapezoidal rule; includes filtering of 
ill-posed initial conditions and automatic selection of 
Jacobi-Nevrton iteration or Newton iteration. 

LSODE 

EPISODE 

VariaUe-step, variable-order backward-differentiation 
method with a generalized Newton iteration*. 

CHEMEQ 

Vai iable-step, second-order predictor-corrector method with 
an asymptotic integration formula for stiff equations. 

DASCRU 

Variable-step, fourth-order, explicit Runge-Kutta-Merson 
solver, | 


* Other options are included in these packages. 


f 

Figure 3 Summary of Methods Studied 
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Two Problems Describing Adiabatic, Constant Pressure Chemical Reactions 


Test Problem 1: 


Combustion of a Mixture of 33% CO and 67% ^ with 100% Theoretical Air 
(taken from reference 11) 

12 Reactions 

11 Species + Temperature 

Pressure = 10 atm. 

Initial Temperature = 1000 K 
Reaction Duration: 1 ms 


Test Problem 2: 


Combustion of a Stoichiometric Mixture of H_ and Air 
(taken from reference 9) 

30 Reactions 

15 Species + Temperature 

Pressure = 2 atm. 

Initial Mixture Temperature = 1500 K 
Reaction Duration: 1 ms 


Both Problems Include All Three Regimes of Combustion: 
Induction, Heat Release and Equilibration. 


Figure 4 Test Problems 
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Figure 5 - Variation with time ol temperature Figure 6 - Variation with time ol temperature 

and species mole fractions lor test problem 1. and species mole tractions tor test problem 2. 

Solution generated with LSOOE-B and EPS • Solution generated with LSOOE-B and EPS • 
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Method A 


For adiabatic, constant-pressure combustion reaction, energy conservation gives 


n h = h = constant (1) 

i i o 


h.j = molal-specific enthalpy of species i (J/kmol) 
h Q = mass-specific enthalpy of mixture (J/kg) 

In this method, equation (1) was solved for the temperature using a Newton- 
Raphson iteration technique. 


NS 

Z 

i=l 


where, 

and 


Method B 


Differentiation of eq. (1) with respect to temperature (T) gives 


dT 

dt 


NS 

.Yi h i 

i^j 

NS 

Z n. c 
1-1 1 P l 


(2) 


where, c p is the constant-pressure specific heat of species i (J/kmol K). 


In this method, the temperature was evaluated by integrating equation (2). 


Figure 7 Evaluation of Temperature 
(for LSODE, EPISODE, and CHEMEQ) 
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HO (s) 

1 

CPU ( s ) 

1 

o 

0.786 

* 

1 

o 

H 

0.703 

10" 7 

0.791 

10-* 

7.91 

10-* 

8.04 

10-*° 

0.772 


Figure 10 Example of Effect of Initial Steplength (HO) 
on Work Required by EPISODE-A (EPS = 10“») 
for Test Problem 2 


An approximate expression for AT — the maximum allowable temperature 
change allowed before the reaction rate constants k. and k . are updated — was 
derived by requiring that the maximum relative errlir in th4 resultant reaction 
rates does not exceed the local relative error tolerance (EPS) required 
of the numerical solution. The approximation for AT is given by 


AT 


EPS-T 


max 

E. 

E . 

7^- + N. 

• ._li + w 

j 

RT j 

’ RT -j 


( 3 ) 


where, T is the current temperature, the bars j denote absolute value, and 
the maximum is taken over all forward and reverse reactions. 


Figure 11 Approximation for AT 
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Method 

Test Problem 
1 2 

GCKP84 

0.85 

1.73 

CREK1D 

0.23 

1.04 

LSODE-ft 

0.31* 

0.52* 

LSODE-B 

0.29* 

0.51* 

EPISODE-ft 

0.75* 

0.54* 

EPISODE-B 

0.70* 

0.67 

CHEHEQ-ft 

6.41* 

13.6* 

ClEltEQ-B 

5.69* 

12.3* 


* method incorporated eq • (3) 


Figure 12 Minimum CPU Time (in seconds on IBM 370/3033 computer) 
Required for the Test Problems 
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